{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sys\n",
    "sys.path.append('..')\n",
    "import mesh, periodic_homogenization, tensors, mode_viewer, numpy as np\n",
    "from tri_mesh_viewer import TriMeshViewer as Viewer"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Load and analyze a 2D microstructure using periodic homogenization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m = mesh.Mesh('../../examples/meshes/2D_microstructure.msh', degree=2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "view = Viewer(m)\n",
    "view.showWireframe()\n",
    "view.avoidRedrawFlicker = True\n",
    "view.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Run periodic homogenization for a base material with E=200, nu=0.35\n",
    "Cbase = tensors.ElasticityTensor2D(200, 0.35)\n",
    "homogenizationResult = periodic_homogenization.homogenize(m, Cbase, orthotropicCell=True)\n",
    "print(homogenizationResult.Ch)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Visualize deformations/stresses of the elastic metamaterial."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compute deformation and strain of homogenized material's minimum energy eigenmode\n",
    "hr = homogenizationResult\n",
    "u, e = periodic_homogenization.probe(m, hr, tensors.SymmetricMatrix(hr.Ch.computeEigenstrains().eigenstrains[0]))\n",
    "stress = Cbase.doubleContract(e)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "view.update(scalarField=stress.vonMises())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "view.update(scalarField=[np.max(s[0]) for s in stress.eigendecomposition()]) # maximum principal stress"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "view.update(scalarField=[np.max(np.abs(s[0])) for s in stress.eigendecomposition()]) # maximum abs principal stress"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "view.update(vectorField=u) # maximum abs principal stress"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Visualize the lowest energy eigenmode with an animation.\n",
    "mview = mode_viewer.ModeViewer(m, u.flatten(), [hr.Ch.computeEigenstrains().eigenvalues[0]], amplitude=0.1)\n",
    "mview.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Load and analyze a 3D microstructure using periodic homogenization.\n",
    "The microstructure we analyze here has reflectional symmetries, and for efficiency, we analyze only the \"orthotropic base cell\" (i.e., the positive octant) using MeshFEM's \"ortho base cell homogenization\" features. At the end of this notebook we demonstrate that this analysis produces the same result as running standard periodic homogenization on the full cell."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m3d = mesh.Mesh('../../examples/meshes/3D_microstructure_orthocell.msh', degree=2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "view3d = Viewer(m3d)\n",
    "view3d.showWireframe(True)\n",
    "view3d.avoidRedrawFlicker = True\n",
    "view3d.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Cbase = tensors.ElasticityTensor3D(200, 0.35)\n",
    "periodic_homogenization.benchmark_reset()\n",
    "hr3d = periodic_homogenization.homogenize(m3d, Cbase, orthotropicCell=True)\n",
    "periodic_homogenization.benchmark_report()\n",
    "print(hr3d.Ch)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compute deformation and strain of homogenized material's minimum energy eigenmode\n",
    "ed = hr3d.Ch.computeEigenstrains()\n",
    "eigenstrains, eigenvalues = ed.eigenstrains, ed.eigenvalues\n",
    "u3d, e3d = periodic_homogenization.probe(m3d, hr3d, tensors.SymmetricMatrix(eigenstrains[0]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Visualize von Mises stress under the minimum energy eigenmode.\n",
    "view3d.wireframeMaterial().color = '#57B'\n",
    "view3d.update(scalarField=Cbase.doubleContract(e3d).vonMises())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Use a different colormap for visualizing stress.\n",
    "import vis, matplotlib\n",
    "sfield = vis.fields.ScalarField(m3d, Cbase.doubleContract(e3d).vonMises(), colormap=matplotlib.cm.viridis)\n",
    "view3d.showWireframe(False)\n",
    "view3d.update(scalarField=sfield)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compute metamaterial's deformation for all its energy eigemodes.\n",
    "modes = np.column_stack([periodic_homogenization.probe(m3d, hr3d, tensors.SymmetricMatrix(es))[0].ravel() for es in eigenstrains])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Visualize one of these modal displacements as a vector field.\n",
    "view3d.update(vectorField=modes[:,3].reshape((-1, 3)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Visualize all the vibrational modes of the 3D elastic metamaterial."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mview3d = mode_viewer.ModeViewer(m3d, modes, eigenvalues, amplitude=0.1)\n",
    "mview3d.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Compare against periodic homogenization on the full unit cell.\n",
    "We demonstrate that the more efficient homogenization routine operating on the positive octant of a structure with reflectional symmetries produces the same result as running traditional periodic homogenization on the full unit cell."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m3d_full = mesh.Mesh('../../examples/meshes/3D_microstructure.msh', degree=2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "periodic_homogenization.benchmark_reset()\n",
    "hr3d_full = periodic_homogenization.homogenize(m3d_full, Cbase)\n",
    "periodic_homogenization.benchmark_report()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(f'Traditional periodic homogenization result:\\n{hr3d_full.Ch}')\n",
    "print(f'Orthotropic base cell homogenization result:\\n{hr3d.Ch}')\n",
    "print(f'Moduli discrepancy: {np.linalg.norm(np.array(hr3d_full.Ch.getOrthotropicParameters()) - np.array(hr3d.Ch.getOrthotropicParameters())) / np.linalg.norm(hr3d.Ch.getOrthotropicParameters())}')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
